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We expand on a previous study of fronts in finite particle number reaction-diffusion systems in the 
presence of a reaction rate gradient in the direction of the front motion. We study the system via 
reaction-diffusion equations, using the expedient of a cutoff in the reaction rate below some critical 
density to capture the essential role of fluctuations in the system. For large density, the velocity is 
large, which allows for an approximate analytic treatment. We derive an analytic approximation 
for the front velocity dependence on bulk particle density, showing that the velocity indeed diverges 
in the infinite density limit. The form in which diffusion is implemented, namely nearest-neighbor 
hopping on a lattice, is seen to have an essential impact on the nature of the divergence. 

I. INTRODUCTION 

The propagation of fronts connecting different macroscopic states is a common occurrence in many non-equilibrium 
systems 0. Familiar examples range from solidification 2] to chemical reaction dynamics such as flames Q and 
to the spatial spread of infections |J| through a susceptible population. Previous work by many authors has shown 
that a useful way to classify such fronts is via the stability properties of the state being invaded. In fact, surprising 
differences, with regard to the selection of the speed of propagation [jj, the rate of approach to that speed the 
sensitivity to finite-particle number fluctuations 0, 0] , and the stability to 2d undulations , exist between fronts 
that propagate into metastable versus linearly unstable states. 

In a previous work [TTj . we addressed the question of the dynamics a new type of front, that which exists in a system 
containing a reaction-rate gradient in the direction of front motion |12|. Our staring point was a simple infection 
model A + B — ► 2 A on a Id lattice (with spacing a) with equal A and B hoppin g ra tes 3 ; this process leads in the 
mean-field limit to a spatially discrete version of the well-known Fisher equation 
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fa = r<t>i(l - (pi) + -y (<pi+i - 20; + . (1) 

Here propagation is into the linearly unstable <j) — state and 4> is just the number of A particles at a site, normalized 
by N. We then assumed that the reaction rate would be a linear function of space, increasing in the direction of 
propagation, This type of gradient would be a natural consequence of spatial inhomogeneity, or could be imposed 
via a temperature gradient in a chemical reaction analog. Also, this type of system arises naturally in models of 
Darwinian evolution [LH Il5j|. (where fitness x is the independent variable; the birth-rate, akin to the reaction-rate 
here, is proportional to fitness). The naive equation describing such a model is the Fisher equation Q with a reaction 
strength r = r a {x) varying linearly in space |16| 

r a (x) =r + ax . (2) 

This model gives rise to an accelerating front. We also introduced a quasi-static version of the model |17| . wherein 
the reaction rate function moves along with the front: 

r q (x) = max(r mi „, f + a(x - x f (t))) , (3) 

with Xf is the instantaneous front position, the precise definition of which we will discuss later. The minimum reaction 
rate r m in is introduced so as to stabilize the bulk = 1 state, and plays no essential role in the following. This quasi- 
static problem should lead to a translation-invariant front with fixed speed v q (fo,a). Although important on its 
own, one might also try to view the quasi-static problem as a zeroth-order approximation to the original model, (the 
absolute gradient case), where by ignoring the acceleration, one obtains an adiabatic approximation to the velocity 
v(t; ro, a) ~ v q (fo(t), a) with ro(t) = ro + axj(t). In both models, the reaction rate gradient was seen to enhance 
tremendously the role of fluctuations, to the extent that the naive treatment via a reaction-diffusion, or mean-field, 
equation gave rise to "irregular" behavior completely at odds with the original stochastic model. In particular, the 
reaction-diffusion system exhibited an extreme sensitivity to initial conditions not present in the stochastic model. 
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Furthermore, the quasi-static version of the reaction-diffusion system exhibited a front which accelerated without end, 
whereas the stochastic version of the model always achieved an asymptotic constant velocity steady state. 

To get some insight into the stochastic model, we employed a heuristic approach in which we mimic the leading-order 
effect of finite population number fluctuations by introducing a cutoff in the mean- field equation (MFE) Q \m\jj%. 
This cutoff replaces r(x) by zero if the density <j> falls below k/N for some 0(1) constant k; this change in the reaction 
term prevents the leading edge from spreading too far, too fast. This idea has proven its reliability in the Fisher 
system with constant reaction rate where it correctly predicts the aforementioned anomalous effects 0- Simulation 
results showed that the cutoff MFE does a quantitatively accurate job of tracking the actual front dynamics. We 
then used the cutoff MFE to study the front velocity, at a fixed spatial position, as a function of N. This was done 
both for the absolute gradient model and for the corresponding quasistatic model. From the data, we concluded that 
both models exhibit velocities which increase, evidently without bound, as a function of N, which is of course radically 
different than what had been encountered in the previous classes of propagating fronts. Thus the cutoff treatment 
succeeded in showing why the long-time dynamics of the stochastic model is not at all correctly described by the naive 
reaction-diffusion system. In addition, the cutoff theory had the physically reasonable property, again in accord with 
the stochastic system, that at small enough TV, the velocity could be approximated by just taking a cutoff version of 
the usual Fisher equation result for a fixed reaction rate rp = r$ + ax, i.e. neglecting the reaction-rate gradient across 
the front. This is so because the effective interfacial width, the distance over which the particle density drops from 
its bulk value 0(1) to its cutoff value 0(1/ N) scales as log N; hence one can neglect the gradient if a log AT is small. 
The naive reaction-diffusion system, however, due to its ever increasing interface width, always feels the reaction rate 
gradient and never in this adiabatic regime. 

Given the highly unusual velocity results, an analytic treatment of the cutoff system at large iV is clearly worthwhile. 
A very telegraphic version of this analysis, as applied to the quasi-static model, was presented in Ref. IllL The purpose 
of this paper is to present this analysis in detail. The continuum problem is treated first in Sec. ^ A treatment of 
the dependence of the velocity on the "base" reaction rate, tq is presented in the next section. In Sec. II VI we redo the 
continuum problem via a WKB treatment, developing the methods which will prove necessary for the lattice problem. 
The lattice problem is attacked in Sec. producing the controlling (geometrical optics) WKB approximation, whose 
properties are then investigated. The full leading order (physical optics) WKB solution is obtained in Sec. IVII This 
is matched to the solution past the cutoff in Sec. IVIII completing the analysis of the model. A summary and some 
concluding remarks then follow. 

II. THE CONTINUUM PROBLEM 

In this paper we study the steady-state motion of fronts in the quasi-static version of our model. On the lattice, 
the solution has the " Slepyan" traveling wave form [2(j : 

(t>i(t) = <p(t-ia/v) , (4) 

where the field <f>i at each lattice site, i, has the same history, shifted in time. The equation of motion then becomes 
a differential-difference equation for <j), which we write as a function of the variable x = —vt : 

= -zMx + a) + cj>(x - a) - 2<f>{x)) + vtf + r{x)(cf> - 2 )6»(0 - 1/N e ) . (5) 

The cutoff in the reaction rate sets is when the density drops below some fraction k of one particle per site, so that 
4> < k/N = 1/N e . We found that k = 0.25 yields excellent quantitative agreement with the stochastic model. In 
this quasi-static model, the reaction rate is also only a function of the comoving variable x: 

r(x) = max(r m „,r + ax) , (6) 

where we have chosen the origin of time such that the front position Xf, defined by </> x , = 1/2, is located at i = at 
t = 0. This definition is simplest one for the deterministic problem posed by the cutoff MFE; other conventions are 
more convenient for simulation studies of the stochastic model > but this merely corresponds to a slight change of 

In the spatial continuum limit, o — > 0, this steady-state becomes a standard differential equation: 

= D<j)" + v(j)' + r{x)(4> - <j) 2 )e(ct) - l/N e ) (7) 

In this section we treat this simpler problem, returning to the lattice version in Sec. Ivl 

We want to solve the problem for large N e , where as we have noted we expect the velocity to be large. If v is indeed 
large, then it appears that the diffusion term is negligible in comparison, and can be dropped. We will see that this 
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is in fact valid as long as x is not too large, including the entire "bulk" region of the solution where <j> is O(l). We 
then get 



v<f>' = -r{x){4> - 4> 2 )0{4> - 1/N e ) 
with the solution satisfying </>(0) — 1/2 given by 

(r x + ax 2 /2)/v 



In 



1 - 



7"min(«£ ^min)/^ H~ (^min 

+ ax min /2)/v 



(8) 



(9) 



where the upper term is valid for x > x m - m , and the lower term is valid for x < x mul . x m { a is the point where the 
minimum reaction rate is reached, 1 + rQax m i n — r m i n . If we assume for the moment that the solution is valid up to 
x c , where the cutoff sets in (i.e. <j){x c ) = 1/N e ), then all we have to do is solve for x > x c and match. The solution 
there is 
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-v{x—Xc) / D 



To do the matching, it is enough to use the small <j> approximation of Eq. JHJ, namely 

i _ ~(r a x+ax 2 /2)/v 



(10) 



(11) 



The matching of <p and </)' /<j) at x c then gives 



1 

v 
D 



-(rox c +ax c /2)/v 



Tq + ax c 



(12) 



two equations for the two unknowns v and x c . For large N e , both of these are large and we obtain the approximate 
solution 



so that 



\nN e = ax 2 J{2v) 
v 2 / D = ax c 



(2D 2 alnN e ) 1/3 
(^) 1/3 (ln^/3 



(13) 



(14) 



We can now check our assumption concerning the irrelevance of diffusion for x < x c . Using Eq. (|llf) . we find that 



D<t>" 



= D 



(ro + ax) 2 jv 2 + a/v 
(r + ax) 



(15) 



For x of order 1, this is of order 1/v and is indeed small. However, x c is large, of order v 2 , so that here the ratio is 
order 1 and diffusion can no longer be ignored. However, since the ratio is of order 1, and not large, the scaling given 
by Eq. ^]is correct, just not the numerical coefficient. 

To incorporate diffusion for x < x c , we can linearize the equation since <fi is already small in this region. We get 



= Dip" + v4> + (r + ax)cj> 
Up to a similarity transformation, this is the Airy equation, with the general solution 



-vx/1D 



where 
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FIG. 1: Exact quasi-static velocity compared to various approximations for the spatial continuum case, a — > 0, with a — 0.3, 
D = 1, rrj = 1. 



We need to match this to the diffusionless solution Eq. (|llfl for 1 <C x <C x c , where the arguments of the Airy 
functions are large and positive. Doing this, we find that B must be set equal to zero, since Bi((j — x)/S) decreases for 
increasing x and so enhances the fast descent of the exponential factor. The Ai term on the other hand increases with 
increasing x and cancels out the fast exponential, leaving the desired slow exponential of the bulk solution. Matching 
to the bulk solution, we find 



We are now again at a position to perform the match at x c - The matching equations are 

u Xc /2D AM [T 



(19) 
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Ai' 



7 — X c 



rJAi( 
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(20) 



Examining the second of this set of equations, we see that the second term on the left must be large, which we can 
arrange if the denominator is small; i.e., Ai is close to its first zero. To leading order in v, 7 ~ v 2 /{AD a) S> 1 and so 
x c w 7. Then, to leading order, we get 



2D 3\5 



2AD 2 a 



so that 



{2AD 2 a\nN e ) 



1/3 



(21) 



(22) 



so that indeed incorporating the diffusion just modified the prcfactor. Plotting together the exact numerical solution, 
obtained from a straightforward Euler initial value integration in time of Eq. , with an appropriately small a and 
r q {x) as the reaction term, with the numerical solution to our analytic matching formula Eq. (|20|l and our asymptotic 
scaling solution Eq. I|22|) . we see that our matching formula agrees extremely well with the exact velocity. Even for 
the large values of In N e considered here, however, the leading order formula is not very impressive. The next order 
term can be calculated and gives 



v = 2{D 2 a) 1/3 [(3 In A^) 1 / 3 + £ (31n Ay~ 1/3 



(23) 



where £0 = —2.3381 is the location of the first zero of the Airy function. Thus, although the correction does decrease 
with N e , it does so extremely slowly. This improved approximation is also presented in the figure, and does extremely 
well. 
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FIG. 2: Comparison of the predicted ro dependence, Eq. I|25jl . and numerical results for the spatial continuum limit a — > 0, 
with a — 0.1, D = 1. Also included is the Fisher velocity, Eq. 126H . where the graphs for the two JV e 's coincide on the scale of 
the figure. 



III. DEPENDENCE ON ro 



To the order considered, ro has not entered into the calculated velocity. We can calculate the leading ro dependence 
by expanding Eq. 121|) to quadratic order in ro: 

lnN e = ^-*(l) 3/2 ~-£--J± (24) 
2D 3\SJ 24D 2 a 2va v ; 



which induces a correction Av to the velocity of 



4D 2 r 2 . . 

Av = ^ (25) 



so that v increasing quadratically with ro to leading order. Note that the shift is small, of order 1/v even for ro of 
order v. Also, ro cannot be taken to be of order v 2 , since then diffusion is relevant in the bulk. In Fig. 0a comparison 
between formula (|25|l and numerical results is shown. The initial quadratic dependence is clearly seen. For large 
enough ro our formula fails. Indeed for very large ro, the effect of the reaction gradient is suppressed, and the velocity 
should approach that of the (cutoff) Fisher equation with rate r 0] 



IV. THE CONTINUUM PROBLEM, A LA WKB 

The fact that we could solve the linear problem exactly obscures the fact that most of the structure of the problem 
comes from the asymptotic properties of the solution. In fact, we can get essentially everything we require via a WKB 
treatment. Writing (f> — e s , we get to leading order 

DS' 2 + vS' + ay = (27) 

where we have written y = x + r$/a and the derivatives are w.r.t. y. Equation (|27fl defines S' implicitly in term of y. 
For small y, we get S' = —ay/v, so that S — So — ay 2 /(2v). Matching to the bulk solution gives So = r o/(^ va ) ~ 0j 
if we take ro to be order 1. Later, we will see what happens when ro is of order v. Now what is critical, as we saw 
above, is the turning point, since beyond this point S' turns complex, 4> starts to oscillate, and so hits zero, which 
allows us to match to the post-cutoff solution. The turning point is given by the discriminant condition, which we 
can write as 



JL[DS' 2 + vS> + ay]=0 



(28) 
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Si = -t^ (29) 
2D v ; 

Solving for the turning point y» gives us y* — v 2 /(4Da) which is consistent with the solution given above in sectionlTTl 
where the turning point occurs where the argument of the Ai function is zero. However, we do not actually need the 
value of the turning point, just that of S 1 there, namely Since, as we verify later, the turning point is close to the 
zero of the solution, the dominant contribution to the value of is e s * . This is given by 



S* = I"* dyS'= I"' dS'S'^- 

rSL 







dS' S'[-(2DS' + v)/a] = (30) 



Note that we did not need an explicit expression for S'(y), which is good, since in the lattice case we won't have such 
an expression. This calculation is already enough to give us the leading asymptotics, since to leading exponential 
order <j>{x c ) = e s ' , or 5* = — rniV e , exactly what we got above. The origin of the correction lies in the fact that the 
zero lies a small distance beyond the turning point, namely £q6, a result we need the Airy equation to derive. The 
real part of S' is fixed at beyond the turning point, so S at the zero is S* — £o<W (2-D), which we need to set equal 
to — In N e . This gives us the correction derived above. 

V. THE LATTICE PROBLEM 

Now we are in a position to return to our lattice problem, Eq. (jSJ). As above, in the bulk diffusion is irrelevant 
and the solution is the same as before. Close to the turning point, we linearize and expand S(x ± a) (even though we 
can't expand <j)(x ± a)), 0,0] and the WKB equation is 

0= ^-smh 2 (aS'/2) + vS / + ay (31) 
or 

Already at this point, we get a nontrivial result. We can de-dimensionalize this equation by introducing T — a /IS, 
z = y/£, £ = v/(aa) so that the equation reads 

0= — sinh 2 (772) +T' + z (32) 
va 

where the derivative is now w.r.t. z. Thus, S (i.e. \nN e ) scales like D/(aa 3 ) times a function of the dimensionlcss 
parameter va/D, so that the results for all a, (for a given k and D) should lie on a universal curve. Furthermore, we 
see that a is a singular perturbation as far as the large velocity limit goes |15). since no matter how small a is, the 
parameter va/ D eventually goes to infinity. 

Returning to Eq. (|31|) . the turning point is given by the discriminant equation 

2D 

0= — sinhfaSi) + v (33) 
a 

which gives 

, 11 v 2 a 2 va\ , 

S * = a 1 »{V 1 + lD 2 --2D) ^ 

Again, we need to calculate the change in S from y = to the turning point y*. This is given as above by 

s - - C is ' s '% 

s '* 2D 

dS' S'[-( smh(aS')+v)/a] 

o a 

2D 2D v(^') 2 

— Si cosh(a^) + smh(aSi) - (35) 
aa z aa 6 2a 
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This then is the leading order WKB answer. Again, to get the correction, we need to examine the vicinity of the 
turning point more closely. We write 



<t>(y) 



= P siv 



(36) 



In the vicinity of the turning point, this removes the large variation of 4> between lattice points, leaving us free to 
Taylor expand the rest. Substituting this into Eq. ©, we get 







D r 

D 

~n2 



ip(y + a) 



-aS' 



fp(y - a) - 2ip(y) + vS'^y) + vij)' (y) + ayi)) 



e aS - (ip + aip' + aV/2) + e~ a5 * (if) - atp' + a 2 ^" /2) - 2ip + vSitp + vip' + ayip 



D cosh(aSl)ip" + a{y — y*)ip 



(37) 



Again, we get an Airy equation. This gives us the distance from the turning point to the zero of ip, which is t;o5 a 
where 



S a = (a/(D coshes;)))- 1 / 3 



(38) 



is the length scale of the Airy equation for the lattice problem. This gives us an additional contribution of — S'^oS a 
to S. Again, the solution for the velocity is just In N e = —S, so that 



ln7V„ 



9D 9D 

cosh«) - sinh«) + v(Si) 2 /2 



DcoshiaS'J 



-1/3 



(39) 



Let us examine the various limits of the this expression. First, the continuum limit, av/D <C 1. Then = —v/2D, 
as in the continuum calculation, so that aS' < 1. Then, 5* = -(2D(S'J 3 /3 + v(Sl) 2 /2)/a = -v 3 / (24D 2 a), also 
exactly as in the continuum calculation. The correction term is — S'^oi^/D)^ 1 ^ 3 , which also agrees. 

Now, as we mentioned above, for any finite a, av/D is eventually large for sufficiently large N e . Then, 
5* = - \n(va/D)/a. This gives 5* = v/(a 2 a)(— \a{va/D) + 1 + ln 2 (wa/D)/2). Now, for very large va/D, 
5* ps —v\n(va/D)/(a 2 a). However, this is only valid for ln(va/D) ^> 2. In fact, it is a reasonable (20%) ap- 
proximation only for \n(va/D) bigger than 10, so that v would be unreasonably large. Thus, a strict asymptotic 
expansion is of no use whatsoever. Over the range 5 < x < 11, an excellent approximation of ln 2 (cc)/2 + 1 — ln(ic) is 
x/7.4 (see Fig. 0. 

Thus, in the relevant range, S*» ps v 2 / (7 AaaD) . Thus, while formally, to leading order 



lniV e re v/(a 2 a)\n 2 (va/D) 

or equivalently 

v w a 2 a\n(N e )/\n 2 (a 3 a\n(N e )/D) 
this is true only for astronomically large N e . More useful, though phenomenological, is 

v ps V '7 AaaD \nN e 



(40) 
(41) 
(42) 
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FIG. 4: Velocity vs. ln(JV e ) for lattice spacing a = 1, a — 0.1, D — ro = 1 from simulation, together with the various analytic 
approximations. The curve labeled "Leading Order" represents Eq. 1351 and that labeled "Correction" represents Eq. 1391 . 
The curve labeled "Effective" represents Eq. 14311 . 
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FIG. 5: Velocity vs. ln(JV e ) for lattice spacings a = 1, 0.5, 0.25, and 0, a — 0.1, D — ro — 1 from simulation, together with the 
analytic approximation Eq. 1391 

Thus, the velocity increases with a, D and N e . Including the correction term, we get the effective approximation 

v re ^7AaaD\nN e - 1.76a 1/3 £> 2/3 (ln(/V e ))" 1/3 ln(7.4a 3 cvln(V e )/L>) (43) 

Fig. 0] presents the case a = 1. We see that for v's bigger than 4, the corrected approximation is excellent. The 
leading order approximation, however, is poor even for ln(iV e ) as unreasonably large as 100. The simplified effective 
approximation Eq. (|43|l is as good as the full corrected approximation for this range of N e . The extremely simple Eq. 
(14211 is as good as the leading order approximation. In Fig. |5j we present data for a number of values of a, ranging 
from to 1, along with our analytic prediction. Again, the agreement is very good. 

In Fig. we show the dependence on a, the gradient of the reaction rate. We see that the rise in v is quite steep 
at first, and then tapers off to an much slower rise. It should be noted how much an effect the correction term has, 
especially at larger a. Nowhere does it look simply proportional to a, as would naively appear from the leading order 
calculation, Eq. 

In Fig. we show the dependence of v on the diffusion constant, D. Again, the rise in v is steep at small D, 
and grows essentially linear for large D. As v/D is a decreasing function of D, for large D the continuum limit is 
eventually valid. We see in fact that the continuum approximation, Eq. Il2.'i[l works quite well over the entire range of 
D. 



VI. NEXT-ORDER WKB 

To go further, we need to both improve our WKB solution and to extend our solution to the region x > x c . Note 
that it is easy to check that we do not need to reconsider the diffusionless solution in the bulk, as the first-order 
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FIG. 6: Velocity vs. a for lattice spacing a = 1, D = ro = 1, ln(jV e ) = 25, from simulation, together with the leading-order 
and corrected analytic approximations, as in Fig. [I] 
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FIG. 7: Velocity vs. D for lattice spacing a = 1, a = 0.1, ro = 1, ln(7V e ) = 25, from simulation, together with the leading-order 
and corrected analytic approximations, as in Fig. [I] Also included is the continuum approximation, Eq. 123H . 



correction to that solution is lower order (in terms of its ultimate effect on the velocity) than the correction we derive 
here. We first consider the next-order WKB solution. We write <j> — e So+Sl , and get the next order WKB equation 



2D 



with the solution 



so that 



aS[ sinh(aS£) + —S(j cosh(a^) 



1 



+ vS[ = 



4> « (av/2D + Bmh.(aS' ))-V*e Cl+So 



(44) 

(45) 
(46) 



We have to match the solution to the bulk solution, which is approximately e~~^r . The matching area is defined by 
the requirement that on the one hand diffusion be irrelevant, so that —5*0 « ay/v <C hxv, or y <^ vhxv, and on the 
other, <j> is small, implying y 2 /v » 1, or y ^> w 1 / 2 . Thus to do the matching we can take y ~ 0(v). Neglecting the 
sinh, we get from Eq. (|46ll : 



and the WKB solution is, 



= e So 



2 \2DJ 



2D 

1 + smh(aS") 

va 



-1/2 



(47) 



(48) 
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We now need to match the WKB solution to the Airy solution of 4> — e s ^ v ^ y *^, where ip = C2^i((y* — y)/S a ). 
First we find the matching region. Clearly we need the argument of the Airy function to be large. Thus gives us that 
!/,-!/> 8 a . In the lattice limit, this reduces to - ?/ > u _1 / 3 , while in the continuum limit, we get !/* - j/ > 1. 
Near the turning point, Sq ~ S 1 ^ + Sa 3 ^ 2 y/y* — y. This is valid as long as aSa 3 ^ 2 \/y* — y -C 1, or y* — y -C S^/a 2 . 
Thus, we can match as long as 8 a <C S^/a 2 , which is uniformly true in the large v limit. Approximating the Airy 
solution, we get 

« C 2 e s '*(y~y^Ai((y*-y)/5 a ) 

« _ y) -mc 2 e s '-(y-y'h-i«v*-y^ 3/2 m 

Approximating the WKB solution, we get 

2D N - 1/2 

va 



2D 
1 + 

va 

2D 



smh(aS' (y))J e So{y) 
sinh(a5^) + a cosh(aS^ )\JSa 3 (y* - y) 



-1/2 



^cosh(aS:)V<5a 3 (y* -y) ) e s '- s '^-y^ 3/2 ^ 3/2 CO) 



Matching these two gives 



C 2 = — / m e s ' (51) 

It is easy to verify that this agrees in the a — > limit with the direct continuum calculation, Eq. H19J1 . 

VII. MATCHING TO x > x c 

Finally, we have to actually match to the solution for x > x c - For x > x c , we clearly cannot expand ip as we have 
done, since the falloff of is much faster than exp(— S' r x). We can understand what happens by considering doing 
the expansion of ip to one more order. The higher derivative induces a boundary layer at x c , which serves to insure 
continuity of the second derivative, but leaves the lower-order derivatives untouched. As we expand to higher and 
higher order, there are more and more boundary-layer modes, which we have to match. The correct way to do the 
matching then is via a Wiener-Hopf (WH) procedure. To do the WH, we consider our problem in the immediate 
vicinity of x c . Here, we can approximate ay by ay c , since as we have seen, y c is large. Now we have a constant 
coefficient difference equation, which we can solve via WH. 

Let us do this first for the continuum problem for practice, since here we know the correct answer. The approximate 
equation is 

D<j>" + v<t>' + ay c e(y c -y)<f> = (52) 
Writing <j) ~ 4>l + 4>r and Fourier transforming, we get 

- Dq 2 (fa +<p R )+ ivq{4> L +</)r)+ ay D <f> L = (53) 

or equivalently 

(-Dq 2 + ivq)<p R = -(-Dq 2 + ivq + ay c )(f> L (54) 

The right-hand operator R(q) = —Dq 2 + ivq has two zeros, one at q — and one at q — iv/D, but 4>R nas a pole 
only at q — qn = iv/D. The left-hand operator L(q) = —Dq 2 + ivq + ay c has two essentially degenerate zeros, close 
to q = iv/(2D), both of which are represented as poles in fa. Thus, we rewrite the equation as follows 

/ s, L(q)(q-q R ) 

(q - qn)fa< = 777 \ — <i>L (55) 

R(q) 
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Now, the left-hand side of the equation has no zeros or poles below the line = v/D, and the right-hand side of the 
equation has no zero or poles above, so they must both be equal to a constant, C . Thus, 



C 



•>R = 



q - qr 

CR(q) 



Cq 



L{q){q-qR) (q - q+){q - q~) 
where q± are the two nearly degenerate roots of L(q), 

o± = — ± A 
q± 2D 



with A = ^ v 2D Day " small. Fourier transforming back, we get 



! >R 



iCe -v(y-y c )/D 



= iC 



<1+ 



iq+(v—Vc) _|_ 



c iq- (v—Vc) 



q+-q- q--q+ 

Examining (f>R, we find that iC = 1/N e . Turning to we found above that y c ~ v 2 /(AaD) — £o<5, so 



A « y/-a^ 6/D 



Now, as long as (y — j/ c )A <C 1, we can write 



-v{y- Vc )/{2D) r 



We now have to match this to the Airy function. Putting y c = v 2 /(AD a) — £q5 — y\, we find that 

-v(y-yc)/2D Ai '(£o)yi ( x _ y-y- 



= Ae~ v ^' 2D t 

Comparing the two results we find y± — 2D/v, and 

Ae -„ W2D 2DAi'(£ ) _ 1 



2/1 



v5 



(56) 



(57) 



(58) 



(59) 



(60) 



(61) 



(62) 



lnAL 



D 



3 (ADa) 3 / 2 V a 
v 3 v£ 



-1/2 



8D 2 a 2D 2 A 



ln(Ai'(Co)2ev / 2^^ 1/3 a 1 /%- 1 / 2 



-In 



/2_DAi^o 
I, vS 



2AD 2 a 2(D 2 a) 1 / 3 

This can be shown to agree with the direct asymptotic solution of Eq. 1201 

Now we do the same for the lattice problem. Again, after setting ay = ay c , the equation has the form 



where 



R{q)<f> R = -L(q)c» L 



r./ \ 4D . 2 aq 

L(q) = R(q) + ay c 



(63) 



(64) 



(65) 



R(q) has two pure imaginary roots, one at q = and the other with positive imaginary part, kr. Since we cannot 
allow <pR to become negative, the dominant solution for large x must be controlled by this positive imaginary root, 
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leading to a pure exponential decay. Only this root and the complex roots which decay faster (Qq > kr) are then 
permissible. As in the continuum, since y c is close to L(q) has a pair of almost degenerate roots, q± with small 
real parts and a positive imaginary part smaller than kr. To match to the At solution, this must be the dominant 
contribution for large negative y — y c , so the only acceptable roots are those with < Qg±. We therefore decompose 
R{q) and L(q) into two factors, one with its zeros below ikr, which we label by a U B" superscript, and the other with 
its zeros above (or equal) , which we label by "[/" . Formally, 

R(q) = R u (q)R B (q) 

L(q) = L U (q)L B (q) (66) 

where 

R U (q) = 
R B (q) = 
L u {q) = 

We have chosen to explicitly break out the factors relating to q± in L B since they will play an essential role in the 
following, and the factor ivq in R B so that the correct behavior at q = is maintained. Then, via the standard 
Wiener-Hopf argument, 

4>R(q) = C—tf 

MQ) = (68) 

It is easiest to proceed if we regularize the problem by effectively discretizing time, replacing the ivq term in the 
operators L, R by ivnt/a sin(qa/n t ) where n t is some large integer. Then, the operators become polynomials of order 
2n t in the variable e lqa ^ nt . Further study reveals that R u then has n t — 1 zeros, R B has n t + 1, L u has n t — 2 zeros, 
and L B has nt + 2. Thus, (f>L behaves as q _1 for large q, and so 

q B 

<Pl(Vc) = -i hm q(j> L (q) = -Cvq+q- TT (69) 

so that 

C = -W^U q -¥ (70) 
Nevq+q- J -- L q£ t 

so that 

4 > L\qj = tt/ w rll b~ (71) 

If y is not too close to y c , only the two dominant modes, which we have labeled q± = —iS'^ ± A survive, and we get 

My) « ^ II ~^~1'V ^> ( 1 + £ f^^B + - y c ) ) (72) 

This is seen to reproduce the continuum results above when a — > 0. Now we must match Eq. (17211 to our Airy function 
solution to Eq. (|37|l . Actually, to the order we are working, we must take into consideration the first lattice correction 
to the Airy equation, namely 

DasmhiaS'J y,, + Dcosh(aS ^ + a(y _ y ^ = _ (73) 
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FIG. 8: Large v analytic approximations for j/i vs. velocity for a — 0.1, ro = 1, D = 1, both for lattice spacing a = 1, and for 
the continuum limit, together with results from simulation. 



The (first order) approximate solution to this is 

m *C 2 (1+ ml (y - y.) j Ai + 



a tanh(aS^) 
6^ 



Substituting y c = y* — £,oS a — Dx, with yi/S a <C 1, we find 



Hy) « c2e s '*^-y-h s '*(y-yJ^M 

5a 



atanh(aS*) 

2/1 h t (2/ - 2/cj 



Thus gives us 



2/i 



atanh(aS^) 
6 



E 



5** 



< 3 



Si 



(74) 



(75) 



(76) 



A graph of y\ as a function of v is presented in Fig. |S| together with the continuum result. We see that whereas y\ 
falls with v in the continuum, due to the lattice correction to the Airy equation, the lattice yi approaches a constant 
for large v. In fact, at large v, all the g's can be calculated analytically, and the some performed. This calculation 
shows that to leading order, both the sum over the right and left modes approaches 1/2, with the difference vanishing 
as v — > oo. The lattice y\ is dominated by the lattice correction to the Airy equation, which approaches the constant 
a 1 6 for large v. Included in this figure is a comparison between analytical and numerical results. We see that as 
expected the analytic result approaches the numerical results are v increases, being quite accurate everywhere. 

The numerical results presented in this graph, in contrast to those presented throughout the rest of this paper, were 
not obtained through direct numerical simulation of the time-dependent equations, due to the high accuracy required 
to perform the comparison with the theory. Our direct numerical simulations were performed via a straightforward 
Euler simulation, which is only first-order accurate in the time step. Extremely small time-steps would have been 
required to obtain the requisite accuracy. Instead, we solved the linearized steady-state equation directly. As opposed 
to v(N e ), which requires a full nonlinear solution, y\(v) is determined solely by the linearized equation. The procedure 
we employed was as follows: The solution past the cutoff was written as a linear superposition of the allowed modes, 
corresponding to the roots g^. The solution to the left of some conveniently chosen yg was written as a linear 
superposition of modes, with the reaction rate set at the constant value r{yi) = r{y c ) — a(y c — ye). The steady-state 
equations between ye and y c were written as a banded matrix, acting on the three sets of unknowns: the coefficients 
of the pre-ye modes; the values of the field between ye and y c ; and the coefficients of the post-j/ c modes. This matrix 
depends on the two parameters v and r{y c ). Now the pre-y^ modes include two real modes, corresponding to the 
two solutions of the quadratic equation for S'. In order to match the Airy function behavior, the faster of these two 
modes must not be present. This is an eigenvalue condition of r(y c ) for a given v. From r(y c ) we can back out j/i, 
as presented in the figure. This procedure converges quadratically in the discretization Ay. The convergence with 
respect to ye is exponentially rapid, and presented no problem. A comparison of the answers obtained in this manner 
with that obtained by direct numerical simulation, at low values of v for which the latter calculation was feasible, 
verified the validity of this alternate approach. This new approach also sheds an interesting light on the selection 
problem inherent in the linearized steady-state equation. 
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16 

FIG. 9: Ratio of N e as predicted by Eq. 11771 to the exact N e from numerical simulation for D = 1, a = 0.1, a — 1. Also shown 
is the ratio of N e as predicted from the lower-order result, Eq. (13911 to the exact N e . 

All that remains is to put everything together and construct the full approximation for the velocity as a function 
of the cutoff, N e . As we have seen, in the matching region, the WH 4>l has the same functional form as the Airy 
solution, once y\ is picked appropriately as described above, the two solutions differing only in normalization. Setting 
the normalization factors equal then fixes 1/N e in terms of our previously calculated C2, (see Eq. (1511) ). Doing this 
gives 

lniV e = -S. + S'M Q 5 a + Vl ) - In (^r 1 ) - In R Zfs'-^ ) ' (7?) 

In Fig. we present the ratio of N e as predicted by this formula to the results of numerical simulation. We see 
that the ratio appears to approach unity as v increases. Together with this is shown the ratio of N e from the lower 
order formula, Eq. Ij39j) ■ This formula, while it does better at small v, is seen to diverge from the exact answer with 
increasing v. 

VIII. CONCLUSIONS 

In summary, we have presented an analytical study of the velocity of Fisher fronts in the presence of a gradient. This 
study exploits the fact that the velocity diverges as the local density of reactants increases. This divergence is one of the 
signposts of the extreme sensitivity to fluctuations of this class of models. One of the most surprising consequences of 
this sensitivity is the different order of the divergence in the continuum versus the lattice model; whereas the velocity 
of the front in the continuum limit diverges as lnV 3 (iV),on the lattice, the velocity effectively diverges as vlnJV. 
The relative insensitivity to the details of the matching to the post-cutoff regime is another characteristic feature of 
this problem. It is clear, for example, that the leading order results are completely independent of the post-cutoff 
dynamics. Even the next order, which formally does depend on matching to the solution past x c , is in fact only very 
weakly modified by the "R" modes arising from this region. This lack of strong dependence is no doubt a major part 
of the reason that the phenomenological cutoff theory works as well as it does in describing the stochastic model. 

Obviously, the cutoff MFE approach cannot capture any of the truly stochastic features of the original Markov 
model. Thus, the next step in our overall program for understanding fronts in gradients must involve adding back in 
the residual effects of finite particle number fluctuations to the cutoff theory. Exactly how to do this is already unclear 
in the simpler case of the Fisher equation front, where it has proven difficult to come up with a simple explanation for 
the numerically determined front diffusion constant. jUI^. The first question to be answered for the gradient case 
is whether the front can be described as simply diffusing (albeit with an anomalous diffusion constant) or whether the 
fluctuation effects perhaps lead to even stronger stochasticity. We hope to report on this issue in a future publication. 
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